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■ Variations of the sunspot number are important indicators of the solar activity cycles. The sunspot formation is a result of a dynamo 
1 process inside the Sun, which is far from being understood. We use simple dynamical models of the dynamo process to simulate the 
magnetic field evolution and investigate general properties of the sunspot number variations during the solar cycles. We have found that 
' the classical Parker's model with a standard kinetic helicity quenching cannot represent the typical profiles of the solar-cycle variations of 
T— H ' the sunspot number, and also does not give chaotic solutions. For modeling of the solar cycle properties we use a nonlinear dynamo model 
of Kleeorin and Ruzmaikin (1982), which takes into account dynamics of the turbulent magnetic helicity. We have obtained a series of 
I i , periodic and chaotic solutions for different layers of the convective zone. The solutions qualitatively reproduce some basic observational 
, features of the solar cycle properties, in particular, the relationship between the growth time and the cycle amplitude. Also, on the longer 
. time scale the dynamo model with the magnetic helicity has intermittent solutions, which may be important for modeling long-term 
I ■ variations of the solar cycles. 

o ■ 

I 1 Introduction 

\ Investigation of solar activity has a long history. The 11-year sunspot cycles were discovered by Schwabe 
1 (1844). Detailed observations of sunspots locations on the solar disk and their sizes provided the latitudinal 
£^ . distribution during the cycle, known as the "butterfly diagram". Discovered by Hale (1919) the flip-flop 
I | phenomenon of the global solar magnetic field with a period of 22 years showed the connection between 
(T) ■ variations of the sunspots number and the magnetic field evolution. Taking into account the Hale's law 
, Bracewell (1953) proposed to include the information about the periodical reversals of the magnetic field in 
O | the sunspot number series by assigning alternating positive and negative signs to the sunspot cycles, thus 
0^ ■ representing the sunspot cycles as a 22-year oscillation. Bracewell suggested that the n onlinear tempora l 
r". \ distribution of the sunspots number can be written in the form of a three-halfs law ( Bracewelll 1 19881 ): 
W(t) = 100 [| Wu n (t) | /83] 3//2 , where Wu n (t) is a new "linearized" sunspot number, which is supposedly 
proportional to a typical strength of the Sun's toroidal magnetic field. 

For an explanation of the magnetic field generation Parker (1955) proposed a simple dynamo model, 
which describes the phenomenon as an action of two factors: the differential rotation and cyclonic convective 
vortices. The mean-field th eory and discovery of the a-effect give us a general description of the process of 
magnetic field generation ( Krause and Radler 1980l ). For recent reviews we refer to Charbonneau (2005) 



and Brandenburg and Subramanian (2005). 

It is known that the dynamo process is ch aracterized by algebraic and dynamic nonlinearities (e.g. 
Kleeorin and Rogachevskii 2007 . Sokoloffll2007l ). The algebraic nonlinear ity can be determined as influence 



of the magnetic field on fluid motions and on the kinetic helicity. This results in quenching of the electro- 
motive force and limits the growth of the magnetic field. The evolution of the small-scale magnetic helicity 
in the turbulent plasma causes a dynamical nonlinearity in the dynamo process. 
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The turbulent helicity conditionally can be divided into two parts: hydrodynamic and magnetic. The 
kinetic helicity describes helical turbulent fluid motions; the magnetic helicity determines the order of 
twisted magnetic field lines. Due to the fact that the kinetic helicity makes the magnetic field small-scaled, 
the back influence on the turbulent fluid motions can restrict the unlimited growth of the magnetic field. In 
the mean-field approach the magnetic helicity is separated into large- and small-scale components. Because 
of the conservation of the total helicity a growth of the large-scale magnetic h elicity due to t he dynamo 
action is compensated by the growth of the small-scale helicity of opposite sign (ISokolofill20o"7n . Thus, the 
small- and large-scale magnetic fields grow together and are mirror-asymmetrical. This means that the 
condition of magnetic helicity conservation is, perhaps, more severe for a restriction of the dynamo action 
than the condition of the energy conservation, which leads to quenching of the kinetic helicity. Since the 
concept of magnetic helicity was introduce into the dynamo theory by Pouquet et al. (1976) many magnetic 
helicity models were suggested. For modeling the solar cycle we choose the formulation of Kleeorin and 
Ruzmaikin (1982), explicitly based on the idea of magnetic helicity conservation. Observational data of 
solar magnetic fields are in a reasonab l e agreement with the idea of the magnetic helicity conservation 
([Kleeorin et a£.ll200l IZhang et a/1120061 . ISokolofj [20071 . ISokoloff et ai1l2007l ). Si milar dynamo models with 



magnetic helicity have been also considered in the context of g alactic dynamos ([Kleeorin et aZ.ll2Q00l ) and 
have included helicity transport by mean flows ( Sur et alll2007l ). which prevents the effect of 'catastrophic' 



helicity quenching (e.g. Brandenburg and Subramanian, 2005). In our model, we don't include the helicity 
transport and consider the helicity quenching with a free parameter. 

In this paper, we consider nonlinear behavior of the Parker's dynamo model (without magnetic helicity) 
and the Kleeorin- Ruzmaikin model (with magnetic helicity), and show that the latter can reproduce the 
qualitative behavior of the sunspot number variations during the solar cycles, including periodic and 
chaotic soluti ons for conditions of the solar connective zone . For simplicity, we use a "low-order model" 
approach (e.g. Weiss et al. 1984 . Sokoloff and Nefyodov 20071 ). reducing the dynamo equations to a simple 
nonlinear dynamical system. In section 2, for consistency, we reproduce the formulations of the Parker's and 
Kleeorin-Ruzmaikin's models, and discuss conditions of linear stability and solar parameters. In section 3, 
we present the numerical solutions and compare these with observed properties of the sunspot number 
variations in solar cycles. 



2 Formulation 



2.1 Parker's migratory dynamo 



In a kinematic approximation the dynamo problem can be described by the induction equation ([Parker 
19551 ) 



<9B 

— = V x (v x B) + 7/ m V 2 B, (1) 

where B is the vector of magnetic field, v is the vector of fluid velocity, r\ m = l/(/j,a) is the molecular 
magnetic diffusivity. The magnetic field, B, and the fluid velocity, v, can be separated in two components 
representing mean and fluctuating (turbulent) parts, or B = (B) + b and v = (v) + u. Here (B) represents 
the averaged over longitude magnetic field, b is the fluctuating part of B, (v) represents mean global-scale 
motions in the Sun, (such as the differential rotation), u is velocity of turbulent convective motions. 
Taking into account that the average of fluctuations is zero ((b) = a nd (u) = 0) fo r the case of isotropic 
turbulence we obtain the following mean-field induction equation (e.g. Moffatt 19781 ) 



^f=Vx «v) x (B) + a (B) - r/V x (B» , (2) 

where rj describes the total magnetic diffusion, which is the sum of the turbulent and molecular magnetic 
diffusivity, rj = rj t + rj m (usually r] m <C rj t ); parameter a is helicity. The first term of the equation describes 
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Figure 1. Local Cartesian reference coordinate system. 

transport of magnetic field lines with fluid, the second term describes the a-effect, and the last term 
determines diffusion and dissipation of the field. 

For describing the average magnetic field, following iParker (jl955h . we choose a local coordinate system, 
xyz, where axis z will represents the radial coordinate, axis y is the azimuthal coordinate and axis x 
coincides with colatitude (figure 1). Effects of sphericity are not included in this model. Hence, the vector 
of the mean field, (B), can be represented as 



(B) = B y (x, y)e y + V x [A(x, y)e y ] 



(3) 



where B y (x,y) is the toroidal component of magnetic field, A(x,y) is the vector-potential of the poloidal 
field. Assuming that (v) = v v (x)e v (rotational component) we can write the dynamical system describing 
Parker's model of the a-dynamo (Parker 19551 ) in the standard form: 



dA 

at 

dB 



aB + r/V 2 A 



at ox 

where G = d {v y ) /dz is the rotational shear. Here for simplicity we omit the subscript of B. 

Assuming that the coefficients are constants and seeking a solution of the model in the form (A, B y ) 
(Aq, Bq) exp[i(fcx — ut)], we find the well-known result that a pure periodic solution exists if 



(4) 
(5) 



D 



aG 
r/ 2 k 3 



(6) 



where D is the so-called "dynamo number". The solutions grow in time for D > 2, and decay for D < 2. As 
described in section 3, in the one-mode approximation the classical Parker's dynamo model even in nonlin- 
ear cases gives only non-chaotic oscillatory solutions, and, therefore, cannot explain the observed variations 
of the sunspot number in the solar cycles. For creating chaotic variations of the magnetic field in the low- 
mode approximation it is necessary to add to the Parker's model a third equat ion describing variations 
of th e magnetic helicity a nd its interaction with the large-scale magnetic field (jKleeorin and Ruzmaikin 
1981 iKleeorin et d\ 1 199.51 ). We note that the absence of chaotic solutions in the Parker's model with a 
simple algebraic quenching may be caused by the low-order t runcation flS chmalz and Stix1ll99lh . and that 
in a general PDE formulation chaotic solutions may appear (jCovas et flllll997l . Il998l i~ 



2.2 The Kleeorin- Ruzmaikin model 



In this se ction for clarity we repro d uce t he derivation of the equation for the magnetic helicity variations 
following Kleeorin and Ruzmaikin ( 19821 ). For this we consider helicity a as a variable consisting of two 
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parts, hydrodynamic (ah) and magnetic (a m ) ( Pouquet et al. 1976 . Kleeorin et al. 19951 ): 



a = au + a m , 

a h = -r (u • (V x u)> /3, a m = r (b • (V x b)> /(12irp), 



(7) 
(8) 



where r ~ Zq/ 7 ?* i s the lifetime of turbulent eddies, Iq is the characteristic scale of turbulent motions in the 
convective zone, and rjt is the turbulent diffusion coefficient. It i s convenient to define the influ ence of the 
magnetic helicity on magnetic field using spectral density x(k) ( Kleeorin and Ruzmaikin 19821 ) 



X= (a-b) = J x(k)dk, 



0) 



where a is the fluctuating part of the vector-potential of magnetic field, A, and A; is a wavevector. 

To derive an equation for the averaged helicity density we multiply the basic induction equation written 
without the differential rotation term 



dB 

~dt 



V x (v x B - r/ m V x B), 



(10) 



by the fluctuating part of the vector potential, a; and also multiply the equation for the vector-potential 

dA 



dt 



v x B — r/ m V x V x A, 



(11) 



by the fluctuating part of magnetic field b, and finally sum these equations. Then, using the averaging we 
obtain 

<A 7 a • ^ + b • ^\ = (a • (V x (v x (B»)> + (a • (V x (v x b))> + (b • (v x <B»> + (12) 



dt 



+ (b • (v x b)) - r/ m «a • (V x (V x b))> + (b • (V x b))» 



Taking into account that b = V x a, after some transformations we obtain the following expression for 
the helicity density 



dx 
dt 



-2([vxb]-(B))-2r ?m (b-Vxb). 



(13) 



Two terms (A [a x [v x (B)]]) and (A[a x [v x b]]), disappear as a result of averaging over the volume. 
Using the m ean-field elect rodynamics approximation and retain only the first two terms for the mean 
electric field (jMoffattlll978h 



e = (v x b) ^ a (B) - <q (V x (B» , 



we obtain 



^ = 2 (p (B) • (V x <B» - a (B) 2 - r/ m (b • V x b)) 



(14) 



(15) 



Then we find a relationship between a m and x- F° r this we assume that in inertial range (ko,k\) the 
last term of (15) can be written as 



(b • (V x b)) = f 1 k\(k)dk, 

J kg 



(16) 
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where the helicity density spectrum, \(k), can be described as (jKleeorin and Ruzmaikinl 11982 
Kleeorin et al\\l993 . IZeldovich et al\\l98$ ): 



X(k) = X- 



1 fko 
k 



9-1 



(17) 



Here k^ 1 is an outer scale of a turbulence, k^ 1 is the scale of the helicity spectrum cutoff. The value of 
parameter q depends on the type of nonlinear interactions in the turbulent plasma. For example, q = 5/3 
describes the case of the Kolmogorov spectrum, and q = 3/2 corresponds to turbulence of interacting 
Alfven waves (Kraichnan's spectrum). In addition, factor q was selected in such a way that coefficient x 
coincides with the mean helicity. Consequently, the magnetic helicity given by (8) becomes 



a, 



12vrp 



(b • (V x b)> 



kl 



hi 



k 2 n(k) X (k)dk, 



(18) 



where r* = 2T$(k/ko) 1 q , T o is relaxation time, averaged over the energy spectrum. 

Taking into account that \x\ = | ( a b) | ~ B 2 jk§ and r\t = (12to/cq) , and substituting (17) into (18) we 
obtain 



a,, 



IX, 



(19) 



where 



1 

18 



u 4-2g 



2-qJ 471/37/4 V 1-1/JfeJ- 



(20) 



and fc* = ki/kQ. 

Finally, multiplying equations (15) and (20) we obtain the expression for variations of the magnetic 
helicity in terms of the mean magnetic field: 



where 



dan. 
dt 



Q_ 

2np 



(B).(Vx(B))-^(B) 2 
V 



(21) 



Q 





(22) 



Equation (21) is written for the case of an uniform turbulent diffusion (rjt = (6toA;q) _1 ), and when the 
magnetic Reynolds number is large, r? ~ rjt. In this model we do not include transport of magnetic field 
by flows. 

For further analysis of the Kleeorin-Ruzmaikin model we transform equations (4)-(5) and (21) in a 
nonlinear dynamical system in non-dimensional variables. Following the approach of Weiss et al. (1984) 
we average the system of equations (4-5 and 21) in a vertical layer to eliminate z-dependence of A and B 
and consider a single Fourier mode propagating in the x-direction assuming A = A(t)e lkx , B = B(t)e lkx ; 
then we get the following system of equations 



— = aB - V k 2 A, — = ikGA - rjk 2 B, 
at at 



dt 



Q_ 

2irp 



-ABk 2 + -(B 2 - k 2 A 2 ) 



(23) 
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This transformation allows us to investigate more easily various nonlinear regimes, from periodic to 
chaotic, and obtain relationships of the basic properties, such as the cycle growth and decay times, duration 
and amplitude. However, we note that the formulation and the interpretation of solutions of the simplified 
system are not straig htforward because it do e s not adequately describes nonlinear coupling of the spatial 
harmonics. Recently, Sokoloff and Nefyodovl ( 2007 ) developed a self-consistent method of reducing the 



dynamo equations, (10) and (11), to a dynamical system. In particular, it follows from their study that 
the dominant modes of the toroidal field in the case of the solar dynamo are described by the harmonics, 
which are antisymmetric with respect to the equator (in accordance with the Hale law), sin(fcx), where x is 
colatitude, and the wavenumber, k, is even: k = 2,4, ... The first k = 2 harmonic has the largest growth rate. 
We retain only this mode in our dynamical model. This nonlinear dynamical system is solved numerically 
by using standard Mathematica software for high-precision integration of potentially stiff problems. 

For interpretation of solutions of the dynamical system in terms of the sunspot number properties we 
use the imaginary part of the toroidal component B(t) because it gives the amplitude of the antisymmetric 
harmonics, and approximate the sunspot number, W, as (ImB) 3 ^ 2 , following Bracewell's suggestion. Of 
course, there might be different definitions of the relationship between characteristics of the magnetic field 
and the sunspot number parameter. We note that the solutions of the dynamical system are qualitatively 
similar for the different harmonics. Nevertheless, we choose the parameters, which correspond to the solar 
situation. 

Making the following substitutions: A = AqA, B = BqB, t = Tot, k = k/r (r is a layer radius), 
Tq = l/{k 2 rf) and a m = aoa m , and taking into account that Aq = B^nk/G we obtain: 



di A - , dB . * - 
— = DB - A, — = iA- B, 
dt dt 



da m 

— j- = -va m + 
dt 



AB — D [B 2 — \A 2 



(24) 



where D = Dqo and a = a n + a m are the non-dimensional dynamo number and helicity, Dq = a^Gr 3 /rj 2 , 
«o = 2Qkv \/G, va is the Alfven speed, v = Tq/T (here we assume that ko is close to k\ then T ~ 
l/(2k 2 r) m ), jkleeorin and EuzmaikirJll982l ^ and A = (k 2 r//G) 2 = Rm~ 2 , k is a characteristic wavelength, 



Rm is the magnetic Reynolds number. In the next section we discuss the ratio of parameters when the 
variations become autonomous. 



2.3 Linear stability 

First, we simplify the third equation of system (24) by neglecting the last term, because A is very small 
for large magnetic Reynolds numbers. Note, that hereafter we omit the hat symbol for all non-dimensional 
variables. Hence (24) becomes 

A 4 = DB-A, A ^ = -iA-B, 
dt dt 

= -"<** + [AB - DB 2 ] . (25) 

Using a linearization procedure, A = A' + A , B = B' + B° and a = a' m + ah, and taking into the account 
that the dynamo number is D = Dqu, we obtain the following equations 

^ = D a h B° -A', *E = L4' - B>, 
dt dt 

= _ va > m + [A°B' + A'B° - 2D a h B'B° - 2D a' m B 2 ] . (26) 
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Then, considering solutions in the form ~ e _ltJ< we obtain the condition of linear stability: |Z?o Q /i| > 2, 
in the approximation of B° = and A = 0. Thus, in the case of l-Doa^j < 2 the magnetic field does 
not grow. For l-Dpa/d = 2 we have a periodical solution like in the Parker's model. The linear theory also 
provides the direction of the migration of the dynamo waves with latitude. These waves travel from higher 
to lower latitudes if D is negative, and in the opposite direction for positive D. Thus, to make the model 
consistent with the butterfly diagram we consider the case of negative D. 



2.4 Solar parameters 

In order to estimate the range of parameters of the Kleeorin-Ruzmaikin model, and for modeling the 
solar cycle we used the standard model of the interior structure of the Sun for the top, bottom and 
middle areas of the convective zone (Table 1). The key parameter of the model is the dynamo number 
D = D^a, because its magnitude determines behavior of the magnetic field. We remind that according to 
the condition of linear stability, Doa n should be greater than 2. Taking into account that Dq = aoGr 3 /r] 2 , 
we assume that r/ ~ rjt, G ~ (v y ) /r, where (v y ) is a typical rotational velocity, r is the radius of convective 
layers. Parameter A determines the influence of vector-potential A on variations of magnetic helicity a m . 
From our estimates it follows that for the solar conditions A < 10 -4 . Consequently, we can neglect the 
term with A. The last non-dimensional parameter, included in equations (24) describes the ratio of 
two characteristic times, To and T. Parameter To = l/(k 2 r]) is estimated using the values of turbulent 
diffusivity r\ t = 1.5H p V conv from the mixing-length model, where H p is the pressure scale height, V conv is 
the convective velocity. In the absence of helicity fluxes the value of the damping parameter v is small. 
However, in reality the helicity fluxes increa se the dissipation rate. This can be modeled as a damping 
term, which increases the effective value of v (Bl ackman and Brandenburg! 20031 ) . Because of this the value 
of v is to some extent uncertain. 



3 Numerical calculations 



3.1 Periodic solutions 



According to the model of Parker (equations 4 and 5), the pure harmonic variations of the magnetic field 
occur for the dynamo-number D = 2. It is also known, that for \D\ > 2 we have solutions with increasing 
amplitude. Following a standa r d procedure, we included a nonlinear it y (alpha-quenching) as ah/ (I + £,B 2 ) 
( Ivanova and Ruzmaikin 1977 . Brandenburg and Subramanianl 200 5), where £ is a quenching parameter, 



which limits the growth of the magnetic field amplitude. Examples of nonlinear solutions for the Parker's 
model are shown in figure 2. However, our numerical calculation showed that in the nonlinear regime the 
variations of the toroidal field (figure 2a) and the sunspot number (figure 2b) are also periodic and similar 
to the classical case of the linear harmonic solution for \D\ = 2. For calculations of the model sunspot 
number we used Bracewell's suggestion (private communication) of the connection between the toroidal 
field strength, B, and the sunspot number, W, in the form of the three-halfs law: W ~ B 3 / 2 . 

We note that the model shows a larger amplitude of variations for the poloidal component than for the 
toroidal field, and that there is a phase shift between them. However, comparing of the observed profiles 

Tabic 1. Parameters of the standard model for the different parts of the solar 
convective zone. 



Parameter 


bottom 


middle 


top* 


radius (r), cm 


5 X 10 10 


6 x 10 10 


7 x 10 10 


density (p), g/cm 3 


2 X lO" 1 


4 x 10~ 2 


2 x 10~ 3 


turbulent diffusivity (rjt), cm 2 /s 


2 X 10 13 


3 x 10 13 


1 x 10 13 


radial velocity gradient (G)**, s^ 1 


5 x 10~ 6 


2 x 10" 6 


-2 x 10~ 6 



* Supergranulation layer. 

** Ca lculated from the helioseismology inversion results of ISchou et al\ 
(1998) 
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a) W b) 




5 10 15 20 25 30 



Figure 2. Solutions (in non-dimensional units) for the toroidal magnetic field (a) and the sunspot number (b) in case the Parker's 
dynamo model for different values of the dynamo number, D, and quenching parameter £ = 10 -3 : D = 3 (solid black curve), D = 5 

(grey curve), D = 7 (dotted curve) and D = 10 (dashed curve). 

of the solar cycles (figure 9) with the solutions for the Parker model (figure 2b), we see that, in general, 
the model solutions have a sinusoidal character and do not describe the basic properties of the sunspot 
profiles with rapid growth and slow decay. 

For the Kleeorin-Ruzmaikin model, given by equations (24), the linear instability condition is also 
\D\ = |a/j-Do| > 2. However, in this case the profile of the periodic solutions is not sinusoidal, and depends 
on the initial conditions, Aq and Bq. For higher initial values the amplitude of the nonlinear oscillations 
in the stationary state is higher. However, the shapes of the oscillation profiles are similar. 

Figure 3 illustrates solutions for the model of Kleeorin-Ruzmaikin, and the corresponding variations of 
the sunspot number for different initial conditions. As mentioned, changes of initial values for magnetic field 
components Aq and Bq leads to very similar profiles. In high amplitude cases, dual peaks may appear in the 
variations of the vector potential, A, of the poloidal field. The evolution of the magnetic helicity represents 
a relatively smooth growth followed by a sharp decay. The helicity has maxima when the toroidal field is 
zero. In these calculations the value of parameter v , which describes damping rate of magnetic helicity and 
depends on the turbulence spectrum and the dissipation though helicity fluxes, is of the order of unity. 
Finally, the variations of the sunspot number, W, with the amplitude increase are characterized by higher 
peaks and shorter rising times (figure 3d). Note that in the sunspot number profile we can recognize the 
well-known general properties of the sunspot number profile with the rapid growth at the beginning of a 
cycle and a slow decrease after the maximum. 

With the increase of la^-Dol (l^/i-Dol > 2) the profile of magnetic field variations continue to deform and 
can become unstable with very steep variations of the magnetic field. The solution can be stable again 
if we enhance the back reaction by increasing the quenching parameter. We use the following quenching 




V/, ... ,m .... wj. .. , sr./. ■ . 
188 190 192 194 196 198 



Figure 3. Variations of the magnetic field for the middle convective zone a^Do = —2: v = 1.28, ct^ = 2.439, Do = —0.82 for different 
initial conditions: Bo = 4i, Ao = — O.Oli (dotted curve), Bo = 4i, Ao = — i (dashed curve) and Bo = 1 + 4i, Aq = — i (black curve): a) 
toroidal component, B; b) vector-potential, A, of the poloidal magnetic field; c) magnetic helicity variations; d) evolution of the model 

sunspot number; e) phase portrait of the magnetic field components. 
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2.5 2.6 2.7 2.8 2.9 

Kinetic helicity, a h p 

Figure 4. Dependence of the amplitude of model sunspot number W on a) kinetic helicity aj, at fixed dynamo number Do 

dynamo number Do at fixed q/, = 2.44. 



-0.82; b) 



formula for the kinetic part of helicity, a^ ( Kleeorin et al. 19951 ) 



a 



1 + tB 2 



+ OLr, 



(27) 



Thus we always have a possibility for selecting £ to obtain periodic nonlinear solutions. 

As mentioned earlier, the change in values of the dynamo number and the kinetic helicity leads to a 
change in the amplitude of the magnetic field variations (figure 4) . With the increase of the kinetic helicity, 
aift, above the linear stability level at a fixed dynamo number (in this example, Do = —0.82) the amplitude 
rapidly grows (figure 4a). However, for a fixed helicity value (e.g. ah = 2.44) and increasing value of 
dynamo number Dq the amplitude rapidly grows and then sharply declines (figure 4b). For other values 
of Dq and ah, satisfying the linear instability condition, the solutions are qualitatively similar. 

Figure 5 shows the relationships between of the amplitude of the sunspot number parameter, the growth 
and decay times and the cycle duration, for different values of a n and Dq. These characteristic times were 
determined from the points of minima and maxima of the model sunspot number, W. The crosses represent 
the amplitude of the periodic field variations for Dq = —0.82 and different values of kinetic helicity, a^. 
The circles correspond to the case of constant ah = 2.44 and different values of dynamo number Do- 
In the first case, the relationship between the cycle amplitude and the growth time is well-defined and 
monotonic. However, in the case of a fixed a^ and varying Dq (circles) the amplitude initially, at small 
Dq, follows the same sequence as in the variable ah case, but then at higher values of \Dq\ (shown by 
bigger circles) the amplitude decreases. The decay time (figure 5b) is longer for higher amplitude cycles. 
The relationship between the amplitude and the cycle duration (figure 5c) is more complex. In the case 
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Figure 5. Relationships between the model sunspot number amplitude, W and a) the cycle growth time, b) the decay time, and c) the 
cycle duration. The circles show a sequence for a fixed value of the kinetic helicity, = 2.44 and the dynamo number varying from -7 
to -0.82. The crosses show the case of fixed Do = —0.82 and varying a^, from 2.44 to 3. The size of the crosses and circles is 
proportional to the corresponding values of Do and a^. Other model parameters correspond to the middle convective zone: 



v = 1.28, A = 1.23 X 10 , and the initial conditions are Aq = — i, Bq 



0. 
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10 Nonlinear dynamical modeling of solar cycles 

of fixed Dq the amplitude and the cycle duration increase with the dynamo number, but initially at lower 
ah, just above the stability value, the amplitude rapidly grows but the duration does not change much. 
With further increase of ah the amplitude grows slowly, but the cycles become longer. In the case of the 
varying dynamo number (circles) the relationship is multivalued, that is the cycles of the same duration 
may have different amplitudes. 



3.2 Chaotic solutions 

The transition of the periodic to chaotic solutions occurs when the dynamo number, |a/i-Do|> increases. 
In the transition regime the cycle amplitude becomes modulated: it slowly increases with time, and then 
suddenly and very sharply declines, and then start growing again (figure 6a). Some examples of the 
transition to the chaotic regime for the sunspot number parameter and the development of chaotic behavior 
are shown in figure 6 for a fixed -Do = —0.82 and increasing values of the kinetic helicity, ah- 

At small deviation from the periodic regime the slow rise of amplitude is followed by sharp recession. 
Nevertheless, each subsequent increase in the amplitude is almost identical to the previous one, or in other 
words, we have in this case, a modulated periodic solution (figures 6 a, b). Also, we note that the shape of 
the profile for the each short cycle is still the same as in the case of the periodic solutions with the rapid 
growth and slow decline. With further increase of the kinetic helicity, the general behavior with growing 
cycle amplitudes and subsequent recessions continues, but the variations of amplitude become increasingly 
more chaotic (figures 6 c, d). The amplitude of the solutions is controlled by the quenching parameter, £, 
which was chosen in figure 6 to obtain the solutions with amplitudes of the same order of magnitude. 

In the case of significant deviations from the condition of linear stability the solutions become chaotic for 
all variables of the dynamical system. Figure 7 shows an example of chaotic variations for the magnetic field 
components, the magnetic helicity and the sunspot number parameter. In the chaotic solutions, the peaks 
of the toroidal magnetic field, B (figure 7a) strongly correlate with the peaks of the vector-potential, A, 
and the magnetic helicity, a m , (figures 7 b, c). The growth of the toroidal field also leads to strengthening 
of the poloidal field and strong fluctuations of the magnetic helicity. The comparison of the variations for 
the toroidal component and vector-potential on the phase portrait (figure 7e) shows behavior resembling 
a strange attractor, when the trajectory converges to a steady state, then deviates from it, and returns 
again. 

The peaks ("cycles") of the chaotic solutions still reproduce the typical profiles of the sunspot number 
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Figure 6. Variations of the model sunspot number, W, illustrating the transition from the periodic to chaotic solutions for conditions 
of the middle convective zone, dynamo number Do = —0.82, v = 1.28, A = 1.23 X 10~ 6 and different values of the kinetic helicity ct^ 
(D a h ): a) a h = 2.58, f = 10~ 3 (-2.116), b) a h = 2.8, £ = 10" 3 (-2.296), c) a h = 3, £ = 1.3 X lO" 3 (-2.46), and d) a h = 3.2, 

£ = 3.9 x 10" 3 (-2.624). 
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Figure 7. Example of the chaotic solution for the middle convective zone (u = 1.28, A = 1.23 X 10 , Do = —0.82, q/, = 3.2, 
£ = 3.9 X 10 — 3 ): a) toroidal field B, b) vector-potential A, c) magnetic helicity ct m , d) model sunspot number W and e) the phase 

portrait of the magnetic components. 

variations described for the periodic solutions. However, the duration and amplitude of the cycles vary. 
For comparison in the figure 8 we show the properties of the different cycles of three chaotic solutions 
for 7_?o = —0.82 and: a n = 2.8 (black circles), = 3 (empty circles) and = 3.2 (stars). In particular, 
we can see that the growth time is shorter for stronger cycle (figure 8a). The decay time, in general, is 
also longer for stronger cycle (figure 8b). However, the relationship between the cycle duration and its 
amplitude is again not very certain (figure 8c). 

It is interesting to note that the similar relationships between the amplitude and the charac- 
teristic times, which we found for both periodic and chaotic solutions of the Kleeorin-Ruzmaikin 
model, occur for the real solar cycles (figures 9 and 10). For comparing the properties of the sim- 
ulated variations of the sunspot number and the solar cycles we used the monthly data of the 
sunspot number (figure 9) for the period of 1750 - 2007 from National Geophysical Data Center 
(ftp://ftp.ngdc.noaa.gov/stp/solar-data/sunspot-numbers/rnonthly .pit) . The growth and decay times were 
measured from the sunspot minimum to the maximum and from the maximum to the next minimum, using 
the same procedure as for the model. The results shown in figure 10 resemble the theoretical relationships 
of figures 5 and 8. In particular, the growth time is shorter for stronger cycles. The decay time and the 
duration do not show a clear correlation with the amplitude. 
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Figure 8. Relationships between the amplitude of the model sunspot number and a) the growth time, b) decay time and c) the cycle 
duration, for Dq = —0.82 and different values of the kinetic helicity: ct^ = 2.8 (black circles), = 3 (empty circles), = 3.2 (stars). 

The time scales are non-dimensional (in units To = l/(k 2 n)). 
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Figure 9. Observed monthly sunspot number series for 1755 - 2007 yrs. from National Geophysical Data Center. 
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Figure 10. Properties of the solar cycles of 1755 - 2007: the relationships between the amplitude of the sunspot number and a) the 
cycle growth time, b) the cycle decay time and c) the cycle duration. The time scales are in years. 



4 Discussion 



We have presented a numerical analysis of simple dynamical models describing the nonlinear behavior of 
two dynamo models, the classical Parker's dynamo model with the standard a-quenching and the model 
of Kleeorin and Ruzmaikin (1982), which takes into account effects of turbulent magnetic helicity. The 
models are analyzed in the low-mode approximation with the goal of representing variations of a sunspot 
number parameter during the solar cycles. We found that the Parker's model does not reproduce the typical 
behavior of the sunspot number with a fast growth and slow decay or obtain a chaotic solution in the 
low-mode approximation, even in strong nonlinear regimes. The analysis of the Kleeorin-Ruzmaikin model 
showed the existence of nonlinear periodic and chaotic solutions for conditions of the solar convective zone. 
For this model we obtained the profiles of the sunspot number variations, which qualitatively reproduce 
the typical profile of the solar cycles. 

It is interesting that the properties of the simulated cycles demonstrate good qualitative agreements with 
the properties of the observed solar cycles. In particular, the relationship between the cycle amplitude and 
the growth time looks similar. Of course, such simplified models cannot pretend to provide a realistic 
description of the solar c ycles. More complicated models such as the models with flux transport (e.g. 
Dikpati and Gilmanl lioOfil ) should be considered for more detailed description of the solar magnetic field 



evolution. 

Nevertheless, these results encourage further theoretical and observational investigations of the dynamo 
model with the turbulent magnetic helicity. We also hope that through the use of a simple dynamical model 
reproducing the basic properties of the sola r cycle it will be come possible to apply the data assimilation 
methods for predicting the solar cycles (e.g. Kitiashvilil 20081 ). 
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